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ABSTRACT 

We present a detailed temporal analysis of a set of hydrodynamic and magnetohydrodynamic (MHD) 
simulations of geometrically-thin {h/r ^ 0.05) black hole accretion disks. The black hole potential is 
approximated by the Paczynski-Wiita pseudo-Newtonian potential. In particular, we use our simu- 
lations to critically assess two widely discussed models for high-frequency quasi-periodic oscillations, 
global oscillation modes (diskoseismology) and parametric resonance instabilities. We find that ini- 
tially disturbed hydrodynamic disks clearly display the trapped global g-mode oscillation predicted 
by linear perturbation theory. In contrast, the sustained turbulence produced in the simulated MHD 
disks by the magneto-rotational instability does not excite these trapped g-modes. We cannot say 
at present whether the MHD turbulence actively damps the hydrodynamic g-mode. Our simulated 
MHD disks also fail to display any indications of a parametric resonance instability between the ver- 
tical and radial epicyclic frequencies. On the other hand, we do see characteristic frequencies at any 
given radius in the disk corresponding to local acoustic waves. We also conduct a blind search for any 
quasi-periodic oscillation in a proxy lightcurve based on the instantaneous mass accretion rate of the 
black hole, and place an upper limit of 2% on the total power in any such feature. We highlight the 
importance of correcting for secular changes in the simulated accretion disk when performing temporal 
analyses. 

Subject headings: accretion: accretion disks — black hole physics — magnetic fields — X-ray: binaries 



1. INTRODUCTION 

Rapid X-ray variability is a ubiquitous characteristic 
of accretion onto black holes. Aperiodic X-ray fiuctua- 
tions are seen from both Galactic Black Hole Binaries 
(GBHBs) and Active Galactic Nuclei (AGN) and, ac- 
counting for the inverse scaling of all relevant frequencies 
with black hole mass, appear to have similar characteris- 
tics (Uttley, McHardy & Vaughan 2005; McHardy et al. 
2006). While it is highly tempting to relate this variabil- 
ity to the magnetohydrodynamic (MHD) turbulence that 
is believed to drive the accretion process (Balbus & Haw- 
ley 1991, 1998), the exact physical processes underlying 
the observed fiuctuations remain mysterious. 

GBHBs also display quasi-periodic oscillations (QPOs) 
in their X-ray lightcurves (see review by McClintock & 
Remillard 2003)^. Of particular interest are the high- 
frequency quasi-periodic oscillations (HFQPOs) that are 
seen in the very-high (or steep power law; McClintock 
& Remillard 2003) state of GBHBs. The HFQPOs have 
quality factors of few-to-10, centroid frequencies of order 
100 Hz and appear to be imprinted on the hard X-ray tail 
of the spectrum rather than the thermal disk emission. 
The fact that their frequencies are stable and at least 
loosely comparable to the orbital frequency at the inner- 
most stable circular orbit (ISCO) around the black hole 
suggests that their properties are set by the relativistic 
portions of the gravitational potential. This gives them 
enormous promise as a diagnostic of black hole mass and 
spin. 

^ Department of Astronomy and the Maryland Astronomy Cen- 
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^ Very recently, the first convincing case of a QPO in an AGN 
was reported by Gierlinski et al. (2008). 



However, the utility of HFQPOs to relativistic astro- 
physics is severely limited by the lack of a compelling the- 
oretical framework in which to interpret measurements 
of the frequencies, quality factors and root-mean-square 
(rms) powers. There exist well defined geodesic fre- 
quencies (i.e., the orbital, radial epicyclic and vertical 
epicyclic frequencies) at any given radius in the accre- 
tion disk. However, these frequencies (as well as all non- 
trivial linear combinations) change with radius, and it 
is not clear why the frequencies of any one particular 
radius would be preferentially displayed in the overall 
power-spectrum. 

HFQPOs are commonly found in pairs with an ap- 
proximate 3:2 frequency ratio, and this has been used to 
suggest that a particular radius is picked out due to a res- 
onance. In the parametric resonance model (Abramow- 
icz & Kluzniak 2001, 2003; Abramowicz et al. 2002, 
2003), there is a resonance in the disk at the radius where 
the radial epicyclic frequency and vertical epicyclic fre- 
quency are in small integer ratios. As discussed below, 
the strongest resonance occurs when these frequencies 
are in 3:2 ratio, at least in the simplest manifestation of 
this model. Other resonance models have been examined 
by RezzoUa et al. (2003a,b) and Kato (2004a,b,c). 

Another interesting possibility is that the HFQPOs are 
global oscillation modes of the accretion disk (i.e., "disko- 
seismic" modes). Global modes have been examined 
analytically (using linear theory) on a hydrodynamic 
background in spacetimes that are pseudo-Newtonian 
(Okazaki, Kato, & Fukue 1987; Nowak & Wagoner 1991, 
1992, 1993; Markovic & Lamb 1998), Schwarzschild 
(Kato & Fukue 1980), and Kerr (Kato 1990, 1991, 1993; 
Kato & Honma 1991; Perez et al. 1997; Silbergleit et 
al. 2001; Wagoner et al. 2001; Ortega-Rodriguez et 
al. 2001). Three classes of mode are recovered corre- 
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spending to pressure modes (p-modes), inertial modes 
(conventionally referred to as g-modes even though the 
restoring force results from rotation or inertia, depend- 
ing on your frame of reference; J.Pringle priv. communi- 
cation) and warping/corrugation modes (c- modes). For 
plausible masses and spins, the fundamental (m = 0) g- 
mode was quickly identified as a good candidate for the 
first HFQPO discovered, the 67 Hz oscillation found in 
the system GRS 1915-^105 (Nowak et al. 1997). Cur- 
rent diskoseismology theory does not provide a natural 
explanation of HFQPO pairs with small-integer ratios; 
however, all present analyses are conducted using linear 
theory whereas these HFQPOs pairs would likely arise 
from mode-coupling that would only be revealed by a 
non-linear analysis. 

Clearly, many open questions concerning the physics 
of X-ray variability remain, including the correct inter- 
pretation of HFQPOs. The current dominant paradigm 
for understanding black hole accretion is that the 
magnetorotational instability (MRI; Balbus & Hawley 
1991) drives powerful MHD turbulence, and correlated 
Maxwell stresses within this turbulence mediate the out- 
ward transport of angular momentum that allows accre- 
tion to proceed. However, the connection between the 
MHD turbulence paradigm and models for the aperi- 
odic and quasi-periodic variability remains highly uncer- 
tain. For example, can the MHD turbulence naturally 
produce the rms-flux relation noted in most black hole 
X-ray lightcurves (Uttley & McHardy 2001) and/or the 
log-normal flux distribution found in Cygnus X-1 (Uttley, 
McHardy & Vaughan 2005)? Arc diskoscismic modes ex- 
cited by turbulent fluctuation (Nowak & Wagoner 1993), 
or does the turbulence act to damp such modes (Arras, 
Blaes & Turner 2006)? Do the Maxwell stresses couple 
radial and vertical motions in such a way as to excite 
parametric resonance instabilities of the type identified 
by Abramowicz & Kluzniak (2001) or any other reso- 
nant phenomena? Does the fact that the HFQPOs are 
imprinted on the high-energy tail provide a fundamen- 
tal clue to their origin, or is it a generic consequence of 
any oscillating thermal accretion disk surrounded by a 
Comptonizing corona (Lehr, Wagoner & Wilms 2000)? 

In this paper, we use a set of global hydrodynamic and 
MHD simulations of geometrically-thin accretion disks in 
a pseudo-Newtonian potential to begin an exploration of 
these issues. Our canonical MHD simulation represents a 
thinner disk, and is run for more orbits, than any previ- 
ously published well-resolved 3-d MHD disk simulation. 
This allows us to conduct a more extensive study of the 
temporal variability of such disks than has previously 
been attempted. In § 2, we give a brief review of the 
theory of both local and global hydrodynamic modes of 
black hole accretion disks, as well as the parametric res- 
onance instability model for HFQPOs. § 3 presents our 
study of ideal (inviscid) hydrodynamic disks, both with 
imposed axisymmetry and in full 3-dimcnsions. We find 
prominent trapped g-modes in the axisymmetric simula- 
tions which remain (albeit with diminished amplitude) 
in the full 3-d case. We then study the MHD case in 
§ 4, where we find that the turbulence excites neither 
the diskoseismic modes nor the parametric resonances 
discussed above. Instead, we find that the turbulence ex- 
cites local hydrodynamic waves of the type elucidated by 
Lubow & Pringle (1993). We discuss our results, includ- 



ing a comparison to previous work, in § 5 and conclude 
in § 6. 

2. THEORETICAL PRELIMINARIES 

Here we provide a brief review of some previously es- 
tablished theoretical results that are pertinent to this 
paper. 

2.1. Local oscillations and waves in accretion disks 

There is a very extensive literature on oscillations and 
waves in accretion disks. Here we focus on just those 
aspects of the field that turn out to be relevant for the 
interpretation of our simulation which have been eluci- 
dated most clearly by Lubow & Pringle (1993). 

Lubow & Pringle (1993; hereafter LP93) studied three 
dimensional wave propagation in accretion disks ignoring 
self-gravity. In the case where one ignores vertical mo- 
tions, they show that radial waves obey the well-known 
dispersion relation 

L0'=K'+Clk', (1) 

where Cs is the sound speed (assumed, in this case, to 
be purely a function of r) and k is the radial epicyclic 
frequency given by = + rdfl^/dr (also see Bin- 
ney & Tremaine 1987). Here, 0(r) is the angular fre- 
quency of the background Keplerian flow. LP93 proceed 
to study the propagation of axisymmetric waves in the 
case where the atmosphere has a locally isothermal ver- 
tical structure. They flnd two types of waves. There are 
low-frequency gravity waves for which 

o<uj <n. (2) 

There are also high-frequency acoustic waves which have 

uj^ > {n-f + l)n^, (3) 

where n = 0, 1, 2, . . . and 7 is the adiabatic index. In the 
special case of purely vertical perturbations, the inequal- 
ity becomes an equality, 

uj"^ ^ {n-f + . (4) 

The n = mode corresponds to a bulk vertical displace- 
ment of the disk and subsequent vertical oscillation at 
the vertical epicyclic frequency which, in the analysis of 
LP93 and in all analyses performed in this paper, coin- 
cides with the orbital frequency. In general, the n — th 
mode has n vertical nodes (i.e., locations where the verti- 
cal velocity perturbation vanishes), and is cither even or 
odd depending on whether n is even or odd, respectively. 

As discussed in § IH our simulations demonstrate the 
effectiveness with which MHD turbulence excites these 
local acoustic waves. 

2.2. Global oscillation modes of an accretion disk 

As discussed in the Introduction, several groups have 
studied the global oscillation of black hole accretion 
disks using linear perturbation theory, identifying three 
classes of normal mode (g-modes, p-modes and c- modes). 
Trapped g-modes have received particular attention as 
a possible source of HFQPO, although the other fam- 
ilies of modes may well be relevant. Here we review 
some of the basic results of these analyses, following 
the approach of Nowak & Wagoner (1991, 1992; here- 
after NW91 and NW92). The NW91 and NW92 analy- 
ses are not fully relativistic, instead employing a pseudo- 
Newtonian potential. Thus, these analyses can be readily 
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compared with our pseudo-Newtonian simulations. We 
also note that full general relativistic MHD simulations 
have typically yielded results for slowly rotating black 
holes that are very similar to those obtained with pseudo- 
Newtonian potentials (e.g., Gammie, McKinney, & Toth 
2003; De Villiers & Hawley 2003). 

NW91 and NW92 use a Lagrangian formalism (Fried- 
man & Schutz 1978) and a WKBJ approximation to de- 
rive the linearized equations describing perturbations of 
an inviscid hydrodynamic thin accretion accretion disk 
about a pure Keplerian background state. They also ini- 
tially examined the special case of purely radial oscilla- 
tions and found the standard dispersion relation of disk 
theory (eqn.[T]). Given that the gravity-modified p-modes 
described by this dispersion relation become evanescent 
when up' < k^, global p- modes can be trapped between 
the ISCO (where k = 0) and the radius at which w — k. 
In practice, however, the "leaky" nature of the ISCO 
would seem to make the trapping of these modes ineffec- 
tive. 

The more general case, including perturbations that 
have vertical as well as radial motions, yields more 
promising results. NW92 examine the linearized equa- 
tions describing the behavior of the scalar potential 

SU EE SP/p, (5) 
where SP is the Eulerian variation in the pressure. They 
show that the linearized equations are approximately 
separable into radial and vertical equations, with the 
separation constant being a slowly varying function of 
r, T(r). The general dispersion relation for these modes 
becomes 

[cj^-jT{r)n^]{u;^~K^)=cj^clk\ (6) 

Assuming that 7X5^^ > k^, non-evanescent solutions ex- 
ist for > jTH.'^ (predominantly radial p-modes) or 
< (predominantly vertical modes). Through 
this analysis, NW92 identify a class of global 5-modes 
that are trapped between two evanescent regions, r < r_ 
and r > r+, where K{r±) — lu. In other words, these 
modes are trapped under the peak of the epicyclic fre- 
quency. They principally focus on the m = (axisym- 
metric) modes, and show that the mode frequency is 
only slightly smaller than the maximum radial epicyclic 
frequency Kmax- Radial harmonics of these modes are 
very closely spaced. Thus, the inner radius at which the 
mode becomes evanescent is still a finite distance (and, 
in plausible settings, several vertical scale heights) from 
the ISCO. This raises the interesting possibility of hav- 
ing appreciable power in such modes without significant 
leakage across the ISCO. 

A major issue, however, is the effect of the turbulent 
MHD background state on these modes. The diskoseis- 
mic mode frequencies are comparable to the frequencies 
characterizing the expected MHD turbulent fluctuations 
(which is very different to the situation in the Sun, for 
example, where the observed helioseismic modes have fre- 
quencies that are four orders of magnitude higher than 
the turbulent turnover frequency). Thus, an MHD tur- 
bulent disk is likely to be a hostile environment for any 
diskoseismic modes. Furthermore, magnetic forces can 
lead to a rather gradual transition in flow properties 
around the ISCO, potentially worsening the leakage of 
the trapped g- modes. Considering these mode destruc- 
tion/suppression mechanisms, it is reasonable to suppose 



that mode survival becomes easier in thinner disks since 
(1) the ratio of the typical turbulent cell size to the ra- 
dial extent of the resonant cavity will decrease with disk 
thickness and (2) there are suggestions that the transi- 
tion in fiow properties around the ISCO is sharper in 
thinner disks (Reynolds & Fabian 2008; Shafee et al. 
2008) thereby producing less mode leakage. This raises 
the possibility that a g-mode can be sustained against 
(or even fed by; Nowak & Wagoner 1993) the turbulence 
in a sufficiently thin disk. 

Previously published global MHD disk simulations 
(e.g., Hawley & Krolik 2001) have modeled fiows as thin 
as h/r ~ 0.1 and have not reported diskoseismic modes. 
However, these authors did not conduct a directed search 
for such modes and hence it is not possible to say that 
the modes were really not present. Careful examination 
of local "shearing-box" simulations have indeed failed to 
find trapped (7-modes associated with MHD turbulence 
(Arras, Blaes & Turner 2006), but this issue has yet to be 
examined in a global thin-disk setting. Searching for and 
characterizing these trapped g-modes in global thin-disk 
simulations will be a major theme of our paper. 

2.3. Parametric resonance 

From the point of view of explaining HFQPOs, the 
need for global oscillation modes is diminished if some 
process does indeed select special radii in the accre- 
tion disk. As discussed in the Introduction, the discov- 
ery of pairs of HFQPOs with small integer ratios has 
prompted several groups to examine resonance models. 
In particular, we shall briefly review the parametric res- 
onance model of Abramowicz & Kluzniak (Abramowicz 
& Kluzniak 2001, 2003; Abramowicz et al. 2002, 2003). 

We begin by considering an accretion disk in which the 
flow deviates only slightly from Keplerian so that the po- 
sition of a fluid element (in spherical polar coordinates) 
is given by 

r{t) =ro + 6r{t), e{t) = + 0(0 - ^t. (7) 

We have specialized to the case of axisymmetric pertur- 
bations. By expanding the resulting equations of motion 
to third order (and wrapping up the non-gravitational 
forces into two unspecified force functions), one finds a 
Mathieu-type equation of motion 

(56*, u + ^'^[l + a cos{Kt)]5e + A (56*, * = 0, (8) 

where a is a small constant that describes the coupling 
between the vertical and radial perturbations, and A is 
a (small) damping constant. Here, we have specialized 
the equations of Abramowicz et al. (2003; hereafter 
A2003) to the case where the vertical epicyclic frequency 
is the same as the orbital frequency. This is valid for a 
non-spinning black hole and, in particular, the pseudo- 
Newtonian potential that we use for the simulations in 
this paper. 

One expects a system described by eqn. [8] to exhibit a 
parametric resonance instability when k, = 2r2/n, where 
n is a non-zero positive integer. Given that black hole 
potentials always have k < fi, the smallest value of n 
for which the resonance condition is obeyed is rt — 3, 
i.e., K — 2il/3. This is expected to be the strongest 
of the set of resonances. It is interesting that the fun- 
damental "test-particle" frequencies at this resonant ra- 
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dius have 3:2 frequency ratio in agreement with observa- 
tions of HFQPO pairs. For the Paczynski-Wiita pseudo- 
Newtonian potential we use in our simulations (Paczyn- 
ski & Wiita 1980, hereafter PW), 
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Hence, the 3:2 resonance occurs at r = 9.2rg. It must be 
noted, however, that a full integration of a toy-problem 
by A2003 reveals that higher-order effects shift the loca- 
tion of the resonance, making the ratio of the epicyclic 
frequencies extremely sensitive to the strength of the cou- 
pling between the r and perturbations. While A2003 
suggest that this sensitivity is a strength of the model, 
allowing application to HFQPO pairs in neutron star 
systems which do not have simple integer ratios, it in- 
evitably diminishes the power of this model to explain 
black hole systems. In addition, there is currently no 
physical model of the coupling between the radial and 
vertical motions. The simulations described in this pa- 
per allow us to assess whether magnetic forces couple 
these motions in such a way as to drive a parametric 
resonance instability. 

3. HYDRODYNAMIC DISKS 

3.1. Initial comments 

In the remainder of this paper, we construct and ana- 
lyze numerical simulations of thin disks in order to exam- 
ine their variability properties, focusing on the presence 
of local and global modes, as well as parametric reso- 
nances. Of course, real accretion disks are believed to 
require at least an MHD-level description in order to cap- 
ture the MRI-driven turbulence that transports angular 
momentum and drives accretion. MHD simulations are 
addressed in the next section. However, it is useful to 
begin with a discussion of ideal hydrodynaniic models in 
order to help isolate the various physical effects present 
in these complex systems. That will be the focus of this 
section. 

In order to allow us to perform a set of simulations 
with modest-to-high resolution, we begin with 2-d (ax- 
isymmetric) hydrodynamic simulations. We expect (and 
indeed show) that these axisymmetric models are well 
suited for studying the fundamental m = g-mode. 
However, once we move to MHD, Cowling's anti-dynamo 
theorem (Cowling 1957) leads to qualitatively different 
behavior in axisymmetric compared with full 3-d sim- 
ulations. Hence, all of the MHD models that we shall 
describe are performed in 3-dimensions. Since we will be 
directly comparing results (e.g., g-vaoAe amplitudes) be- 
tween the hydrodynamic and MHD simulations, we also 
need to perform a "bridging" 3-d hydrodynamic simula- 
tion. 



3.2. Basic simulation set-up 

To make the problem tractable despite the severe reso- 
lution requirements imposed by the geometrical thinness 
of the accretion disk, we choose to focus on only the most 
essential aspects of the physics. From the discussion in 
§ 2, it is clear that the essential aspect of the relativis- 
tic potential that must be captured is the nature of the 
radial epicyclic frequency (e.g., that it goes to zero at 
some finite radius and hence produces an ISCO at that 
radius). In this sense, the PW potential (eqn. [9l) is a 
good approximation for the gravitational field of a non- 
rotating black hole; its ISCO (at 6rg) and marginally 
bound orbit (at 4rg) are both at the same radius as in 
the Schwarzschild geometry. We also simplify the sim- 
ulations by neglecting all radiation processes (radiative 
heating, radiative cooling, radiative transfer, and the dy- 
namical effects of radiation pressure). In place of a full 
energy equation, the gas is given an adiabatic equation 
of state with 7 = 5/3. 

All simulations are performed in cylindrical polar co- 
ordinates (r, z^cj)). The initial condition consists of a disk 
with a constant mid-plane density (po = 1) for r > risco- 
The initial density scale-height of the disk is assumed to 
be constant with radius, implying a sound speed which 
falls off with radius as approximately r~^/^. There are 
two motivations for choosing to model a "constant-h" 
disk; (1) such a choice is well suited to the cylindri- 
cal symmetry of our coordinate grid and (2) according 
to the standard model of Shakura & Sunyaev (1973), 
the radiation-pressure dominated disks of real accreting 
black holes are likely to maintain an approximately con- 
stant scale-height in their innermost regions. In more 
detail, the vertical structure of the disk is given by 

p(r,z) = /3o exp ^-^^ , (12) 



p(r, z) - 



GMhl 



■p{r,z), 



(13) 



{R-2rgYR 

where r is the cylindrical radius, z is the vertical height 
above the disk midplane and R — + z"^. This 
corresponds to an isothermal atmosphere which, when 
hi = h2, is in vertical hydrostatic equilibrium in the PW 
potential. In order to give the disk an initial vertical 
kick, we set hi = I.2/12 (the values of /12 for all of our 
runs are detailed in Table 1). Thus, the initial disk tem- 
perature is ^ 20% too cold for the density and pressure 
run, leading to a vertical collapse and bounce of the disk. 
The initial density is set to zero for r < risco- The initial 
velocity field is everywhere set to 



VGMr 
r - 2r„ ' 



Vr =Vz = 0, 



(14) 



corresponding to rotation on cylinders and pure Keple- 
rian motion for material on the mid-plane. We impose 
zero-gradient outflow boundary conditions on both the 
radial and vertical boundaries of the simulation, i.e., the 
fluid quantities in the ghost zones are set to the values 
of the neighboring active zone, and a "diode" condition 
is imposed on the component of the velocity perpendic- 
ular to the boundary which allows outflow but disallows 
inflow. 

Table 1 details our hydrodynamic simulations. We per- 
form four simulations in which strict axisymmetry is im- 
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Fig. 1. — Snapshots of evolution of the canonical axisymmotric hydrodynamic simulation (HD2d_l) at t = O.STisco (top-left), t = l.OTjgco 
(top-right), t = 5Tisco (bottom-left) and t = 50Tisco (bottom-right), where Tisco is the orbital period at the ISCO. Both the color table 
and contours show the logarithmic density structure of the disk cross-section, with 10 contours per decade of density. A range of densities 
spanning three orders of magnitude are shown. The curved line to the left of each frame represents the event horizon. 



posed at all times {d/d(f> = 0). These axisymmetric sim- 
ulations were performed using the serial ZEUS-3D MHD 
code (Stone & Norman 1992a, b) in its pure hydrody- 
namic axisymmetric mode. In our canonical axisymmet- 
ric run (HD2d_l), we set /i2 = 0.3rg (corresponding to 
h2/r = 0.05 at the ISCO). To model the dynamics of 
the disk in a robust manner, our vertical domain must 
cover many scale-heights; in our canonical run, the verti- 
cal domain is z g (— 5/i2, +5/12). We place 128 uniformly 
spaced grid cells in this vertical direction, giving 13 cells 
per scale-height. This allows us to resolve hydrodynamic 
waves with wavelengths of O.5/12 or greater. The verti- 
cal resolution requirements force us to consider a limited 
radial domain. In the canonical simulation, the radial 
domain extends from 4rg (i.e., well within the plunge re- 
gion) to 16rg and should correspond to the range of radii 
where trapped diskoseismic modes or A2003-type reso- 
nances occur. Tolerating a grid-cell aspect ratio of ~ 2, 
we place 256 uniformly spaced radial cells in this radial 
domain. The simulation was evolved for a time 200risco 
where 

ri,co« 61.6 GM/c3, (15) 
is the orbital period at the ISCO. 



In order to test the robustness of the results discussed 
below to resolution and the limited radial domain, we 
performed a second simulation (IID2d_lhr) in which the 
spatial resolution was doubled (i.e., the size of each voxel 
was halved in both the radial and vertical dimensions) 
and the radial domain extended out to 28rg. We also 
performed two additional axisymmetric simulations cm- 
ploying the same set-up as the canonical simulation but 
with disks that are half (HD2d.2) and twice (HD2d_3) 
the thickness (including appropriate modifications to the 
vertical domain and resolution; see Table 1). 

As discussed above, we perform a 3-dimensional hy- 
drodynamic simulation in order to aid the later inter- 
pretation of the 3-d MHD simulations. The set-up 
of our 3-d run (HD3d_l) is identical to that for the 
canonical 2-d run except that the computational domain 
has a 0-dimension. To reduce computational expense 
while capturing the essential physics, we simulate only a 
A(/) — 30° wedge of the disk using 32 uniformly spaced 
grid cells, imposing periodic boundary conditions on the 
0-boundaries. This 3-dimensional simulation was per- 
formed using an MPI-parallelized version of ZEUS kindly 
supplied to us by Eve Ostriker (and similar to the ZEUS- 
MP code of Vernaleo & Reynolds 2006). 
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TABLE 1 

Summary of the hydrodynamic (HD) and MHD simulations discussed in this paper. Column (1) gives the designation of the 
SIMULATION. Column (2) lists the dimensionality of the simulation. Column (3) gives the fractional disk thickness at the 

ISCO; THE ASTERISK (*) DENOTES THAT THE SIMULATION WAS STARTED WITH AN INITIAL VERTICAL STRUCTURE THAT IS SLIGHTLY OUT OF 
EQUILIBRIUM IN ORDER TO SEED SUBSEQUENT OSCILLATION MODES (aS DESCRIBED IN THE TEXT). COLUMNS (4), (5) AND (6) LIST THE r-, Z- 
AND 0- DOMAINS OF THE SIMULATION BOX. COLUMN (7) GIVES THE NUMBER OF COMPUTATIONAL ZONES WITHIN THE DOMAIN. COLUMN (8) 

LISTS THE RUN TIME OF THE SIMULATED DISK. 



3.3. Axisymmetric hydrodynamic models and the 
recovery of trapped g-modes 

We now discuss the evolution of the axisymmetric hy- 
drodynamic simulations, beginning with the canonical 
simulation. Starting from the initial condition, the disk 
undergoes dynamical timescale variability because of 
pressure gradients which push matter inside of the ISCO. 
The strong transients close to the ISCO launch outward- 
radially directed waves into the disk which break rapidly 
to become rolls. This behavior can be seen in Fig. [TJ 
These high amplitude transients are short lived, how- 
ever, lasting only ~ lOTjsco- At subsequent times, the 
disk settles into a stationary state apart from small am- 
plitude (and decaying) internal oscillations and a very 
weak accretion stream driven by the numerical viscosity. 

In order to study the decay of the hydrodynamic fluc- 
tuations in a more quantitative manner, we start by com- 
puting the quantity, 

K = I pvldV. (16) 
Jt> 

This quantity is particularly well suited to diagnose ver- 
tical oscillations of the disk, and will vanish once the 
hydrodynamic configuration has established a stationary 
state. In order to diagnose the state of the body of the 
disk (i.e., side-stepping issues of the outer radial bound- 
ary or the plunge region), we choose to compute this in- 
tegral over a restricted domain 2? consisting of the annu- 
lus r £ {7rg, 14rg). The time-dependence of K/ maux{K) 
for the canonical axisymmetric simulation (HD2d_l) and 
its high-resolution counterpart (HD2d_lhr) is shown in 
Fig. [2] The rapid initial decline of K{t) is very similar 
for these two simulations and corresponds to the strong 
transients described above. At long times (after about 
t ~ 4 X 10^ GM jc^ K. 70 Tisco) the behavior of these sim- 
ulations begins to deviate. IID2d_l continues to decay 
in an approximately exponential manner K(t) oc e"*/"^", 
where tq « 4 x Vd'^GMjc?. Superposed on this decay 
is a distinct oscillation. This corresponds to (twice) the 
frequency of the trapped g-mode that we shall discuss 
below. The high-resolution version of this simulation 
IID2d_lhr shows very similar behavior except that the 
exponential decay time is longer, tq ~ 6 x Wi^GMjc?. 
This suggests that the decay of these small perturba- 
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Fig. 2. — Change of the quantity K = fj^pv^dV with time. 
The integration domain © is the annulus r G (7, 14). The bottom 
(black) hne is for the canonical-2d run {HD2d_l), lower-middle 
(red) line is for the high-resolution 2-d run (HD2d_lhr), the upper- 
middle (blue) line is for the canonical 3-d run (HD3d_l), and the 
top (green) line is for the canonical MHD run (MHD_1). The dis- 
tinct "ringing" seen in runs HD2d_l an d H D2d_lhr corresponds to 
the axisymmetric g-mode (see Section 13.31 1. In run HD3d_l, non- 
axisymmetric aperiodic structures mask the underlying g-mode 
(see Section ITit . The fact that K{t) for the run MHD.l lies an 
order of magnitude above that for HD3d_l demonstrates that the 
MRI-driven turbulence completely overwhelms the hydrodynamic 
disturbance s see n in the analogous 3-d hydrodynamic simulation 
(see Section l4.2| l. 

tions is due to numerical dissipation which scales approx- 
imately as the square root of the size of the simulation 
grid cells. 

We now study the spatio-temporal variability of the 
disk and, in particular, seek the diskoseismic modes pre- 
dicted by linear theory. Figure [3] shows the mid-plane 
value of Vr on the (r, i)-plane, i.e., the value of the func- 
tion Vr{r,z = 0;t), for run HD2d_l. Note that, outside 
of the ISCO, the average value is {vr) ^ 0.001c so that 
this Figure essentially plots the fluctuation of from 
its mean value. At early times, we see strong wave-like 
disturbances which are generated in the inner parts of 
the disk (at r « 8rg) and propagate to both large and 
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small radii. Although the outer radial boundary con- 
dition is zero-gradient outflow, impedance mismatching 
results in some reflection of these initial strong waves. 
After these initial transients have died out, it can be seen 
that the highest amplitude fluctuations are limited to a 
rather narrow range of radii in the approximate range 
r = 7 ~ 9rg. The fact that these perturbations are essen- 
tially vertical on the {t, r)-plane indicates that they are 
coherent across this radial range, as would be expected 
if we are seeing a global mode. 

This temporal variability can be explored in more de- 
tail using the power spectral density (PSD), deflned as 

P{v) — a\f{iy)\'^ , where a is some normalization constant 
and /{v) is the Fourier transform of the time-sequence 
f{t) under consideration. 



fH = fit)e 



dt. 



(17) 



Note that, for ease of interpretation, most of the PSDs 
presented in this paper will be in terms of the usual 
frequency, i/, rather than angular frequency, lu. Fig- 
ure |4] shows the PSD of the mid-plane pressure as a 
function of r and frequency for HD2d_l and HD2d_lhr. 
This is computed using approximately the final half 
{At = 102.4Tisco ~ 6308GAf/c3) of the simulation in 
order to avoid the initial strong transients. These PSDs 
differ in the range r = 12 — IGr^; run HD2d_lhr has sig- 
nificantly less low frequency noise than run HD2d_l in 
this radial range. We attribute this to effects related to 
the outer radial boundary which is at r = 16rg in run 
HD2d_l but substantially further out (r = 28rg) in run 
HD2d.lhr. 

Inside of r = 12^3, however, the PSDs of these two sim- 
ulations are very similar and we can trust that neither 
the resolution nor the outer radial boundary affect the re- 
sults significantly. We note that simulations in which the 
initial radial density profile of the disk is truncated be- 
fore reaching the outer boundary also produce essentially 
identical results, again giving us confidence that noise in- 
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Fig. 3. — Radial component of velocity Vr on the midplane of 
the disk as a function of radius and time for the the canonical ax- 
isymmetric hydrodynamic simulation (HD2d_l). The linear color 
table extends from radial velocities of Vr = —0.001c (black) to 
Vr = 4-0. 001c (white). Note that, outside of the ISCO, the average 
value is \vr\ <C 0.001c so that this diagram essentially plots the 
fluctuation of Vr from its mean value. 



filtration from the outer boundary is not an important 
issue. In addition to very low frequency noise, the most 
prominent feature of the inner-disk PSD is a vertical 
ridge of enhanced power at « (4 — 5) x 10^'^ c" /GM ex- 



tending from ; 



6.5rg out to r : 



9.5r„. The one dimen- 



sional cuts through this ridge in Fig. ?? demonstrate it to 
have all of the expected properties of a trapped g-mode. 
Firstly, it exists in a rather narrow range of frequencies, 
(4 — 5) X 10"'^ (? jGM, just below the maximum 



radial epicyclic frequency (kj, 



5.52 X 10 



-3 ^3 



IGM). 



Secondly, it is spatially bounded by the radii at which 
the mode frequency becomes equal to the radial epicyclic 
frequency. There is, however, some leaking of the mode 
down to the ISCO. 

We can visualize the eigenmode by producing maps 
of pressure deviations that have been "period-folded" on 
the period corresponding to the peak power in this mode. 
More precisely, we use the last half of the simulation 
to produce maps of the difference between the instan- 
taneous pressure and the time-averaged pressure with a 
sampling rate of 0.2risco- We then sort these maps into 
16 phase bins (based on the period corresponding to the 
peak power in this mode) and average together all maps 
within a given phase bin. The result is shown in Fig. [Sj 
In addition to the g-mode itself, these maps reveal a strik- 
ing "chevron" pattern at large radii. Time-sequences of 
maps reveal these features to be outward-radially travel- 
ing acoustic waves driven by the global g-mode, refracted 
into the upper layers of the disk atmosphere as they prop- 
agate. 

As a final demonstration that we have properly identi- 
fied trapped g-modes in our axisymmetric hydrodynamic 
simulations, we use our set of runs to study the depen- 
dence of the mode frequency on the sound speed in the 
disk. Figure [7] shows the dependence of the mode fre- 
quency on midplane sound speed of the disk (and hence 
disk thickness) at r = 8rg. As expected from analytic 
theory (see Appendix A), the difference between the 
mode frequency and the maximum epicyclic frequency 
depends linearly on the sound speed. There is, however, 
a discrepancy in the slope of this linear relationship ob- 
tained by the simulations and expected from the analytic 
theory. Note that the analytic treatment assumes that 
the gas is strictly isothermal across the whole region of 
interest, a condition that is clearly violated in the simu- 
lated disk. Given the sensitivity of the mode to the ra- 
dial structure of the disk (as demonstrated by the factor 
of 2.5 difference in the analytic results between the two 
pseudo-Newtonian potentials examined in Appendix A), 
the non-isothermality of the gas in the simulated disk can 
readily shift the mode frequency away from the analytic 
value. 

3.4. Extension to 3- dimensional hydrodynamic models 

We begin our analysis of the 3-dimensional hydro- 
dynamic simulations (IID3d_l) by examining the de- 
cay of hydrodynamic perturbations in our canonical 3- 
d run (HD3d_l) using the quantity K{t) as defined in 
eqn.flBl The K{t) behavior for HD3d_l is somewhat dif- 
ferent from the axisymmetric simulations at late times, 
reaching a quasi-steady state at K/ ma.x{K) ~ 3 x 10~^ 
with aperiodic fluctuations rather than continuing a ring- 
ing exponential decay (Fig. [J). Numerical dissipation 
must be just as effective in HD3d_l as compared with 
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Fig. 4. — Midplane PSD for pressure fluctuations in run HD2d_l (left panel) and its high-resolution counterpart HD2d_lhr (right panel). 
Note that we show the radial range r S (4, 16) which is coincident with the full computational domain of HD2d_l but only the inner half of 
the domain for HD2d_lhr (whose full radial domain extends from Arg to 28rg). Also shown are the radial epicyclic frequency (solid white 
line), orbital frequency (dashed) and the n = 1, 2,3 pure vertical p-modes (from left to right dot-dashed lines). The absolute scaling of the 
PSD, as indicated by the color-bar, is arbitrary. 




Frequency (c/GM) r (GM/c) 

Fig. 5. — Left panel : Midplane pressure PSD for run HD2d_l summed up over a range of radii Ar = 0.5rg centered on r = 8rg. The 
dashed line shows the maximum radial epicyclic frequency. Right panel : Integral of the midplane pressure PSD of those frequency bins 
that exceed 5 X 10"^* in the left panel, as a function of radius. Vertical dashed lines show the range of radii for which the mean frequency 
of this peak is less than the radial epicyclic frequency. 



HD2d_l, hence these perturbations must be driven by 
some instabihty. While it is thought that free Keple- 
rian accretion disks are stable to linear hydrodynamic 
perturbations (see Balbus & Hawley 1998; Hantao et al. 
2006), the presence of the simulation boundaries can in- 
troduce true instabilities (associated with reflection from 
the boundaries) as well as uncontrolled numerical noise, 
and these might explain these low level sustained fluctua- 
tions. Visual inspection of density slices in the (r, z) and 
(r, (j)) planes also suggests that non-axisymmetric wave- 
like perturbations interacting with the boundaries are 
responsible for perturbing the 3-d hydrodynamic disk. 
However, a detailed study of the sustained fluctuations 
in the 3-d hydrodynamic disks is beyond the scope of this 
paper. Since these perturbations exist at a low level (par- 



ticularly compared with the MHD turbulence discussed 
in § 4; see Fig.[5]for a comparison of the sustained fluctu- 
ations in the 3-d hydrodynamic run compared with the 
canonical MHD run) and are likely to be driven by our 
boundaries (hence, are not of astrophysical relevance), 
they are not of importance for the principal focus of our 
study. 

The principal issue to be addressed here is the impact 
of the transition to 3-dimensions on the presence of the 
trapped m = g-modes in the simulated disks. We 
might expect the power in these axisymmetric modes 
to be reduced in the 3-d case as the free energy is 
shared with non-axisymmetric modes. This is indeed 
the case. Figure |S] shows the PSD at r = 8rg of the 
last At = 102.4risco of the canonical 3-d run HD3d_l. 
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Fig. 6. — Maps of period-folded pressure deviation (i.e., the difference between the instantaneous pressure and the time-averaged pressure) 
for run HD2d_lhr. Only the last 102.4 orbits of data have been included in order to avoid the transient behavior associated with the initial 
conditions. The folding period corresponds to the peak of the PSD see in Fig. ??. Phases of 0.0 (top-left), 0.25 (top-right), 0.5 (bottom-left) 
and 0.75 (bottom-right) are shown. This, in essence, gives us a direct view of the eigenmode. 

A region of enhanced power is clearly seen in the correct 
range of frequencies v « (4 — 5) x 10~^ c^/GM to be iden- 
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Fig. 7. — Frequency of the g-modes as a function of mid-plane 
sound speed at r = 8rg for runs HD2d_2, HD2d_l, and HD2d_3 
(from left to right). The vertical bars indicate the range of fre- 
quencies over which enhanced power is seen. 
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Fig. 8. — Midplane PSD of azimuthally-averaged pressure for run 
HD3d_l summed up over a range of radii Ar = 0.5rg centered on 
r = 8rg . The dashed line shows the maximum epicyclic frequency. 
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tified with the axisymmetric g-modes studied in § 13.31 
Furthermore, examination of the radially-resolved PSD 
shows that this region of enhanced power is bounded 
by the radial epicyclic frequency in precisely the man- 
ner expected for trapped g-modes. A comparison of the 
absolute values of the PSD across the mode does reveal, 
however, that the mode contains almost an order of mag- 
nitude less power than in the axisymmetric case. 

We note the existence of a narrow but large amplitude 
spike just above a frequency of 7 x 10^^ (? jGM in the 
r — 8rg PSD of this run. The identification of this feature 
is not clear; it is not at the frequency of any expected 
global g- or p-mode. It is, however, confined to a single 
frequency bin, contains only a small amount of power 
and only shows up over a narrow range of radii. It seems 
likely that this is a noise spike. 

4. MAGNETOHYDRODYNAMIC DISKS 

Having gained an understanding of the thin hydrody- 
namic disks, we move onto the more astrophysically rel- 
evant case of MHD disks. In this section, we construct 
MHD simulations of thin accretion disks in a PW poten- 
tial. We then examine the properties of the broad-band 
noise and search for modes in the resulting turbulent 
MHD disks. 

4.1. Simulation set-up 

We simulate geometrically-thin MHD accretion disks 
by building upon our hydrodynamic computational set- 
up described in § 13.21 As discussed in § 13.11 we only 
consider 3-dimensional MHD simulations. 

Table 1 details our set of MHD simulations. Most of 
our discussion will center around run MHD_1 (which we 
shall refer to as our canonical MHD run). In this run, the 
hydrodynamic variables are set-up as for the canonical 
hydrodynamic run (see § 13. 2p except that we begin the 
disk in vertical hydrostatic equilibrium {hi ~ h2 = O.Srg 
corresponding to hi/r = /12/r = 0.05 at the ISCO). An 
initially weak magnetic field is introduced in the form 
of poloidal field loops specified in terms of their vector 
potential A = {Ar, A^, Ac/,) in order to ensure that the 
initial field is divergence free. We choose the explicit 
form for the vector potential, 

A^ = Ao f{r, z) sin , A, = A, ^ 0,(18) 

where Aq is a normalization constant and f{r,z) is an 
envelope function that is unity in the body of the disk 
and then smoothly goes to zero at r = ^isco, r — ''out 
and at a location three pressure scale heights away from 
the midplane of the disk. The use of f{r, z) keeps the 
initial field configuration well away from either the ra- 
dial boundaries of the initial disk configuration or the 
vertical boundaries of the calculation domain. The fi- 
nal multiplicative term produces field reversals with a 
radial wavelength of 5h. This results in a number of 
distinct poloidal loops throughout the disk. The nor- 
malization constant Aq is set to give an initial ratio of 
gas-to-magnetic pressure oi (3 = 10^ in the inner disk. In 
our canonical MHD simulation, this initial condition is 
evolved for a duration of 630Tisco (38800 GM/c^) using 
the MPI-parallelized version of ZEUS described above. 

We supplement the ideal MHD algorithms of ZEUS in 
two ways. Firstly, it is necessary to impose a floor to 



the density field of 10^^ times the initial maximum den- 
sity in order to prevent the numerical integration from 
producing negative densities. This essentially amounts 
to a subtle distributed mass source. The density only 
reaches this floor close to the z-boundary (i.e., many scale 
heights above and below the disk plane). Secondly, we 
implement the prescription of Miller & Stone (2000) to 
include some effects of the displacement current, princi- 
pally forcing the Alfven speed to correctly limit to the 
speed of light as the magnetic fields becomes strong. We 
note that this "Alfven speed limiter" only plays a role 
within small patches of the tenuous magnetized atmo- 
sphere that forms at large vertical distances above and 
below the disk; it never plays a significant role in the 
body of the accretion disk. 

Periodic boundary conditions were imposed on the cj)- 
boundaries, and zero-gradient outflow boundary condi- 
tions were imposed at both the inner and outer radial 
boundaries (Stone & Norman 1992a, b). However, the 
choice of the z-boundary condition for this kind of sim- 
ulation is notoriously problematic. The most physically 
motivated choice would be a free outflow boundary. How- 
ever, as described in Stone et al. (1996), field-hne "snap- 
ping" at these free boundaries can halt such a simulation. 
Indeed, our own test simulations employing zero-gradient 
outflow boundary conditions on the z-boundaries were 
found to be subject to these difficulties, as well as occa- 
sional numerical instabilities appearing to result from an 
interplay of the imposition of the density floor and the 
free boundary. Furthermore, these tests showed that the 
tenuous matter high above the disk mid-plane generally 
flows slowly across these boundaries at very sub-sonic 
and sub-Alfvenic speeds; strictly, this invalidates the use 
of such boundary conditions anyways (since the flow on 
the other side of the boundary should be able to act back 
on the simulation domain). 

We adopt the solution of Stone et al. (1996) and 
choose to employ periodic boundary conditions in the 
z-directions. While this is obviously unphysical in the 
sense that matter cannot leave the simulation domain in 
the vertical direction, it does guarantee mathematically 
reasonable behavior at the boundary (eliminating numer- 
ical instabilities) and, more importantly, appears to have 
no effect on the dynamics of the accretion disk itself (as 
diagnosed through comparisons with our vertical-outflow 
test runs). In order to further isolate the simulated disk 
from the vertical boundaries, we expand the vertical do- 
main (compared to the canonical hydrodynamic simula- 
tion) to z e (-3,3) (i.e., ±10/11). 

We perform four additional simulations aimed at 
demonstrating the robustness of the canonical simula- 
tion. In run MHD_2, we restrict the vertical domain back 
to z e (—1.5,1.5), allowing us to gauge the importance 
of the location of the z-boundaries. Run MHD_2hr has 
an identical set-up to MHD_2 except that the radial and 
vertical resolution is doubled compared with the canon- 
ical run (i.e., the voxel size is reduced by a factor of two 
in each of the radial and vertical directions), and the 
azimuthal domain is doubled to A(f> = 60° at fixed res- 
olution. Due to the factor of 8 increase in number of 
computational cells (and the decrease in the timestep), 
this run was only integrated for a duration of 85risco 
(5236 GM /c^). While the run time is insufficient to con- 
duct a detailed temporal study, a comparison with run 
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MHD_2 does allow us to investigate the effect of both the 
radial/ vertical resolution and the extent of the ^-domain 
on the turbulent state (see below). As discussed below, 
we find that the properties of the simulated disks are 
very similar these two runs. 

Run MHD_3 is similar to MHD_2 except that the outer 
radius is pushed to r = 28rg (at fixed resolution), dou- 
bling the size of the radial domain. Comparing MHD_2 
and MHD_3, we find no evidence that the outer radial 
boundary is affecting the inner disk (r < 12rg) in any 
way. 

4.2. Basic evolution of the MHD disks 

We now discuss the evolution and general properties 
of these MHD disks, centering our discussion around 
run MHD_1. At very early times {t < STisco) strong 
hydrodynamic transients dominate the evolution as ra- 
dial pressure forces drive mass into the region within 
the ISCO. Similarly to the hydrodynamic cases, these 
transients launch outwardly directed axisymmetric waves 
that break to form rolls. These strong hydrodynamic 
transients largely damp away in all but the outermost 
parts of the disk within lOTisco- Concurrently, the MRI 
amplifies the initial magnetic field until the (domain 
wide) volume-averaged ratio of gas-to-magnetic pressure 
peaks at 5 (at t « lOTlsco), at which point most of 

the flow has become turbulent. The entire flow (outside 
of the ISCO) becomes turbulent by t « 20Tisco- Figure[9] 
displays the density field across a vertical slice through 
the accretion disk at various times. 

Between t = lOTisco and t = 20risco, the total mag- 
netic energy in the computational domain declines un- 
til ^ 20. After this, a quasi-steady state seems 
to be achieved where the magnetic field generation by 
the MRI-driven MHD turbulence is balanced by the re- 
moval of magnetic field energy due to numerical recon- 
nection and magnetic buoyancy. The fact that buoyancy 
is playing an important role is revealed by examining 
time-sequences of the strengths of B-field components in 
(r, z) cross-sections of the disk. One clearly sees highly 
magnetized structures being generated close to the mid- 
plane of the disk which then propagate vertically away 
from the midplane^. From this time until the end of the 
simulation a,t t = 630Tisco, MHD turbulence and hence 
accretion is sustained. Over the course of the simulation, 
the disk loses a little more than 60% of its mass (Fig. [T0|) . 

The total magnetic energy and thermal energy undergo 
a slow decline as mass is drained out of the simulated 
disk. During this decline, the volume-averaged plasma-/3 
parameter within the domain remains in the range (/3) ~ 
20 — 30. However, as expected, the P parameter within 
the high-density body of the disk is appreciably higher, 
reaching values of /3 « 70—100. While this is significantly 
larger than the /3 found in global simulations of thicker 
disks (e.g., Hawley & Krohk 2001, 2002), it is in hue with 
what might be expected for thin accretion disks with 
zero net field as diagnosed through local shearing-box 

^ In principle, one could address the importance of magnetic 
buoyancy in removing magnetic energy from the mid-plane regions 
of the disk by comparing the vertical Poynting flux with (numeri- 
cal) reconnection losses. However, energy losses due to numerical 
reconnection are cannot be tracked in our simulation (due to the 
non-conservative nature of ZEUS algorithm) and, hence, a rigorous 
study of this issue is not possible. 



simulations both with and without vertical stratification 
(Stone et al. 1996; Hawley, Gammie, & Balbus 1996; 
Miller & Stone 2000). 

A generic concern in this class of simulation is the affect 
of the z-boundaries. Once it achieves its quasi-steady 
state, our canonical MHD simulation displays a 2-3 order 
of magnitude drop in magnetic pressure, and a 5 order 
of magnitude drop in gas pressure, between the disk and 
the z-boundary. Thus, the disk boundary seems to be 
well isolated from the boundary. Further confidence is 
gained from an examination of run MHD_2 in which the 
z-boundaries have been brought in from z = ±3 to z = 
±1.5. Despite the fact that the magnetic pressure now 
only drops by 1 order of magnitude from the disk to the 
boundary, all of the results from the canonical MHD run 
discussed in this paper are reproduced by MHD_2. More 
precisely, (1) the PSD of the fluid variables f ii4.3.2p are 
essentially indistinguishable, failing to show any evidence 
for global modes but displaying prominent local p-modes, 
(2) the PSD of the mass accretion rate has a broken 
power-law form with indices that differ from those found 
in MHD_1 by A R:! 0.1 (comparable to the la error bars) 
and break frequencies that differ by A(log ti^rcak) ~ 0.05 
(again, comparable to the la error bars). 

Another generic concern with thin disk simulations is 
whether the resolution in the vertical direction is ade- 
quate. The canonical MHD simulation (MHD_1) has 
approximately 13 grid cells spanning a vertical range 
Az = h2, implying that we cannot follow any modes 
with a wavelength smaller than A ~ /12/2. This is just 
sufficient to follow the fastest growing MRI mode (with 
wavelength A 2Trh/P^/'^) if /? < 100. To ensure that 
we have, in fact, achieved adequate resolution, we com- 
pare runs MHD_2 and its high resolution counterpart, 
run MHD_2hr. For example. Fig. [TT] compares the time- 
dependence of the thermal and magnetic energies, as 
well as the height dependence of the plasma-/? parame- 
ter. While there is some indication of increased buoyancy 
driven escape of magnetic fields from the high-resolution 
simulation (as is apparent from the higher value of P at 
intermediate heights), the two simulations generally com- 
pare very well. Thus, we conclude that we have achieved 
adequate numerical resolution. 

4.3. Temporal power spectra of the basic fluid variables 

4.3.1. The importance of correcting for secular changes 
during the simulation 

The length of our canonical MHD run makes it well 
suited for a detailed study of temporal variability. In 
particular, the long stream of simulation data facilitates 
the construction of PSDs. However, as we now discuss, 
significant complications arise in the analysis of these 
MHD simulation as compared with the pure hydrody- 
namic simulations. 

In the case of the hydrodynamic simulations, the back- 
ground (i.e., unperturbed) flow achieves almost a station- 
ary state once the large transients caused by the initial 
conditions have died out. In particular, the lack of angu- 
lar momentum transport within the hydrodynamic disk 
(other than that due to the very small numerical vis- 
cosity) allows the disk to achieve a non-accreting state. 
Density and pressure fluctuations about this background 
state can then be studied via the straightforward con- 
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Fig. 9. — Snapshots of evolution of the canonical MHD simulation (MHD_1) at t = (top-left), t = ITisco (top-right), t = lOTiaco 
(bottom-left) and t = lOOTisco (bottom-right), where Tigco is the orbital period at the ISCO. Both the color table and contours show the 
logarithmic density structure of the disk cross-section, with 10 contours per decade of density. A range of densities spanning three orders 
of magnitude are shown. The curved line to the left of each frame represents the event horizon. 



struction of the PSD. The constant background level does 
not contribute to the power spectrum and hence the PSD 
faithfully characterizes the fluctuations of interest. 

However, even once the initial transients have been dis- 
sipated, MHD disks never achieve this kind of stationary 
background flow. MHD turbulence leads to continued 
accretion that depletes mass from the simulated disks. 
This decline in total mass leads to secular changes (with 
an approximately exponential form) in the density and 
pressure of the background flow. Unless corrected for, 
even a rather slow exponential decay can have a signif- 
icant influence on the PSD of the pressure or density 
fluctuations, severely affecting attempts to characterize 
the properties of the astrophysically relevant fluctuations 
(i.e., the fluctuations that would be present in the ideal 
case of a steady-state disk in which the mass was replaced 
from a large reservoir). 

To see this, consider some quantity whose time-series 
f{t) we extract from the simulation (for example, this 
could be the mid-plane density or pressure at some given 



radius). Let us assume that this can be decomposed as 

f{t) = eit)[l+git)], (19) 

where g{t) is the time-series of the astrophysically- 
interesting fluctuations about some mean state (i.e., 
{g) = 0), and e{t) is a decay function that describes the 
secular change in the background state due to the drain- 
ing of mass from the simulation. It must be noted that 
the decomposition given in eq. [TO] is not completely gen- 
eral; this assumes that the amplitude of the "real" fluc- 
tuations is proportional to the mean background value, 
and that the properties of the fluctuations otherwise re- 
main invariant as the background state decays. This de- 
composition would be valid if the decay of the simulated 
disk simply amounted to a gradual decline in the den- 
sity scale of the simulation while the temperature and 
velocity structures remained unchanged. We shall refer 
to such (simulated) disks as density-invariant disks. This 
does indeed appear to describe our simulated MHD disks 
(i.e., there is little or no corresponding secular change in 
disk thickness or characteristic turbulent velocities). In 
this case, the mass accretion rate will be proportional to 
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Fig. 10. — Normalized total mass as a function of time for run 
MHD.1. 

the density and the decay will hence have an exponential 
form, e = Co e"*/*". 

In the uninteresting case where there are no fluctua- 
tions {g{t) = Oyt), the observed signal is just f{t) = e{t), 
leading to a Fourier Transform (FT) and PSD given by 

^i'^) — ;, — — , and Pe((^) — ,, 7T, 

{l/tn) + iu;' (l/to)2+tj2' 

(20) 

where Ai and A2 are uninteresting normalization con- 
stants. Thus, for u! ^ 1/^0, the PSD of the exponential 
decay goes as Pe{oj) cx uj^^. 

In the more interesting case of non-zero fluctuations, 
the FT of the observed signal is, 

f{uj) = ?(w) + J e{Lj')g{uj ~ uj') duj', (21) 

i.e., the sum of the FT of the exponential decay with the 
convolution of the FTs of the decay and the interesting 
fluctuations. When the fluctuations are small compared 
with the (decaying) background state, as is the case for 
the density and pressure, one can see that the PSD of 
the observed signal will be dominated by the l/w^ as- 
sociated with the decay. Even in the case where the 
fluctuations are large compared with the decay, the PSD 
will still be influenced by the exponential decay due to 
the convolution term in eqn l21l In particular, regions of 
the "real" PSD which are steeper than (including 
the high-frequency wing of any QPO or regions above a 
high-frequency break) will tend to get filled. 

Clearly, we must correct for this decay of the back- 
ground state, and be cognizant of manifestations of any 
remaining uncorrected effects of this decay. This pro- 
cedure plays the same role as the "pre-whitening" em- 
ployed by Schnittman, Krolik, & Hawley (2006; hereafter 
SKH). We proceed by dividing the observed time series 
by a "best-fitting" exponential decay function. Here, the 
time-constant of the exponential decay io is estimated 
via two methods. Firstly, we can set to by the require- 
ment that the starting and final values of the observed 
series are equal [f{t = 0) = f{t = T), where t = T 
corresponds to the end of the time-series]. This is the 



procedure adopted by Schnittman et al. (2006) except 
that they employ a linear form for the decay function. 
Secondly, we can set io via a least-square fit of an expo- 
nential form to the observed time-series. These methods 
give very similar values of and similar corrected PSDs 
when applied to density or pressure fluctuations, as we 
shall now see. 

4.3.2. Power spectra of fluid variables in disk mid-plane 

We now examine the PSD of the azimuthally averaged 
fluid variables (velocities, pressure and density) for the 
simulated MHD disks, as we did with the hydrodynamic 
disks in §0 The top panels of Fig.fTHshow the PSD for 
the radial and vertical components of velocity in the mid- 
plane of the disk, u™"^(r, i) and w™"^(r, i) for a duration 
lasting /Si = 409.6Tisco starting ai t — lOOTisco (i.e., wefl 
after all of the initial transients have dissipated and a 
quasi-steady turbulent state has been established). Note 
that, outside of the ISCO, these velocity components are 
themselves first-order fluctuating quantities (i.e., the ra- 
dial and vertical velocity of the background state is ap- 
proximately zero) and, within the density-invariant disk 
assumption, will have characteristic fluctuation ampli- 
tudes that remain constant throughout the simulation. 
Hence, we do not need to correct these PSDs for the 
effect of the density decay in the simulation. 

It is readily seen that the PSD of the midplane radial 
velocity w""'^ is dominated by the radial epicyclic fre- 
quency from r ^ 7 — Sr^ out to the outer radial bound- 
ary. Inside of r « 7 — 8rg, the PSD shows broad band 
power associated with the transition from turbulent to 
plunging flow. 

The PSD of the midplane vertical velocity w™"^ is quite 
different and can be interpreted in the light of the lo- 
cal fluid oscillations discussed in § 12.11 Below the ra- 
dial epicyclic frequency, there is a broad spectrum of g- 
modes. Above the orbital frequency, the distinct modes 
described by eqn. 3] become apparent; the n = mode 
[uo > fl) and n = 2 mode [uj > 2.0817) can be picked out 
as distinct tracks, and there are hints of higher frequency 
modes as well. The fact that only even-n modes are seen 
is readily understood given that the odd-n modes have 
Vz-nodes at the midplane, as is evident from the power 
deficit centered on the n — I frequency in the Vz PSD. 

The bottom panels of Fig. [T^] show the correspond- 
ing PSD for the density and pressure from the canonical 
simulation. The displayed PSDs have been corrected for 
the decay of the background state by dividing through 
by an exponential curve that best-fits (in a least square 
sense) the domain-averaged density or pressure. Very 
similar results are obtained using an exponential deter- 
mined from just the end points. The radial banding 
of the low-frequency noise is largely the effect of the 
pre-whitening procedure; the full (non-corrected) low- 
frequency power is somewhat greater and the bands cor- 
respond to the residuals between the real decay and the 
simple exponential model. In addition to this low fre- 
quency noise (which is of secondary interest to the in- 
vestigation presented in this paper), both PSDs show an 
enhancement corresponding to the n — 1 vertical p-mode 
(with to > 1.63i7). As expected, there is no enhancement 
in the density or pressure PSDs corresponding to the 
n = or n = 2 p-modes; these have pressure and density 
nodes at the midplane. 
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Fig. 11. — Left panel : Total (domain integrated) magnetic and thermal energies for run MHD_2 (black solid line) and its high-resolution 
counterpart (MHD_2hr; red dashed line). Right panel : Azimuthally-averaged plasma-/? parameter (i.e., the ratio of the thermal to magnetic 
pressure) as a function of vertical height in the disk at r = Sr^. To obtain this plot, data from t = 30 — SSTjaco have been averaged together. 
The solid (black) line shows run MHD_2 whereas the dashed (red) line shows its high-resolution counterpart (MHD_2hr). 



In stark contrast to the hydrodynamic simulation, the 
(decay-corrected) midplane pressure PSD does not show 
the K-bounded vertical "ridge" on the frequency-radius 
plane that is characteristic of trapped g- modes. The ab- 
sence of an excited trapped g-mode is also confirmed by 
examining the PSD at r = 8rg (Fig. I13i normalized such 
that a direct comparison can be made with the 3-d hy- 
drodynamic results of Fig. [5]). It is important to note, 
however, that a trapped g-mode of the strength seen in 
our 3-d hydrodynamic simulation would not stand out 
from the background turbulence. Thus, while it is appar- 
ent that the fundamental trapped g-mode is not strongly 
excited by the turbulence (in the way that the local p- 
modes are, for example), we cannot say whether the g- 
mode is actively damped by the turbulence. 

Finally, we note that there is no indication that the 
parametric resonance instability of Abramowicz & Kluz- 
niak is at work, at least in its simplest form. As discussed 
in § 2, this model predicts its strongest resonance at the 
location where the radial and vertical epicyclic frequen- 
cies are in a 2:3 ratio; this occurs at r = 9.2rg in the PW 
potential. As shown in Fig. [131 the PSD of the pressure 
fluctuations at r = 9.2rg shows no structure associated 
with the local epicyclic frequencies. We have verified 
that a similar conclusion holds true for the PSD of the 
other variables. While it is possible that non- linear ef- 
fects have shifted the location of the resonance inwards 
from r = 9.2rg (A2003), we note that the PSDs of Fig. [12] 
show no obvious radius at which the power at the radial 
epicyclic and the orbital frequencies appears to be locally 
enhanced. 

4.4. Temporal properties of the instantaneous black hole 
accretion rate 

The previous section addressed the temporal proper- 
ties of the fundamental fluid variables through the body 
of the disk. Of course, real observations of accretion disks 
measure the electromagnetic radiation generated by the 
accretion flow. While our simulations miss all of the 
physics relevant for making a first-principles prediction 
of the observed lightcurve, it is still instructive to ana- 



lyze a simple scalar quantity that can be generated from 
the simulation and may be related to the observed hard 
X-ray radiation (since it is the hard X-rays that carry 
the HFQPO signal). 

Here we consider one such proxy for the observed 
lightcurve, the instantaneous mass accretion rate into the 
black hole, 



M= / r{~Vr)pdS, 



(22) 



where the integral is performed over the surface dTZi 
defining the inner radial boundary of the computational 
domain. Figure [Ml shows the raw lightcurve (i.e., not 
corrected for the exponential decay of the disk). In the 
remainder of this section we address the temporal prop- 
erties of this lightcurve. Of particular interest is the pres- 
ence of breaks or QPOs in the PSD of this light curve. 
Hence, we need a quantitative approach by which the 
significance of such features in the PSD can be assessed. 
We begin by discussing our general statistical approach, 
which differs from the one advocated by SKH. 

4.4.1. Analysis method 

Our approach to PSD fitting is predicated on the as- 
sumption (also made by SKH) that at a given frequency 
the power density has an exponential probability distri- 
bution; if the mean is poj then the probability of measur- 
ing a power between p and p -\- dp is 



P(p)dp = — e 
Po 



-P/Po 



dp 



(23) 



Suppose that we have a model that predicts a power 
density Pmod.i for frequency bin i, and that in our MHD 
simulation we actually observe a power density Pobs,i in 
that bin. The likelihood of the data given the model in 
that bin is then 



C^ = (l/Pmod,i) exp(-pobs,i/Pn 



(24) 



The likelihood of the whole power density spectrum given 
the model is the product of the individual likelihoods, but 
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Fig. 12. — PSD for azimuthally-averaged midplane fluid quantities from in run MHD_1. Panels show PSDs for radial velocity (top left), 
vertical velocity (top right), density (bottom left) and pressure (bottom right). The density and pressure variables have been divided by 
the least-square best-fitting exponential function to correct for the secular decay of the simulation, as discussed in § 14.3.11 Also shown are 
the radial epicyclic frequency (solid), orbital frequency (dashed) and the n = 1,2,3 pure vertical p-modes (from left to right dot-dashed 
lines). The absolute scaling of the PSDs, as indicated by the color- bar, is arbitrary. 

described in §4.3.1 in order to determine and then divide 
out the best-fitting exponential decay. 



the log likelihood is typically more useful: 

In/^ = ^ [- lnPmod,i - Pobs,i/Pmod,: 



(25) 



This is the figure of merit for a given model. It can 
therefore be used for standard tasks such as parameter 
estimation and model comparison (e.g., determining if a 
QPO or break is required by the data). 

Note that to maximize the information content, the 
bin size should be the smallest possible, in this case the 
frequency resolution of the raw PSD. As a result, this 
method does not require rebinning to coarser resolution. 
Our approach therefore yields an accurate evaluation of 
one or more precisely specified models, as opposed to the 
broader but less sensitive method of trying to detect a 
signal in a model-independent way. 

We apply a Markov Chain Monte Carlo method to 
search for best fits and establish confidence regions. 
As discussed in § 14.3.11 we correct for secular changes. 
Given the large amplitude fluctuations in the M and S'tot 
curves, the end-point method discussed in § 14. 3. H is not 
appropriate. Hence, we employ the least-square method 



4.4.2. Results 

The mass accretion rate M has sufficient intrinsic vari- 
ability that secular changes in the disk properties have 
only minor effects on the PSD. This is evident from 
Fig. [151 which compares the PSD and Fourier phase of 
the raw M curve as a function of frequency (left panel) 
with these after multiplication by the exponential in time 
that minimizes the overall rms amplitude (right panel). 
The lack of significant differences suggests that inferences 
drawn from the PSD are robust. The comparison shown 
in Fig. [Tni is for the second quarter of run MHD_1 (ap- 
proximately t = 150risco to t = 310risco); analyses of the 
third and fourth quarters also reveal a lack of sensitivity 
to the secular decay of the disk. 

Having established that the secular decay of the disk 
is unimportant for the PSD of M, we will work directly 
with the raw time series from run MHD_1 , without mul- 
tiplication by an exponential function. We analyze sep- 
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Fig. 13. — PSD of the midplane (decay-corrected) pressure for 
Ar = O.Srg wide zones centered on r = 7rg (top, red curve), r = 
8rg (middle, black curve), and r = 9.2rg (bottom, green curve). 
The thin vertical dashed line marks the maximum radial epicyclic 
frequency. A line-segment with a slope of —2 is shown for reference. 
Also shown (heavy vertical blue lines) are the radial and vertical 
epicyclic frequencies at r = 9.2rg where the simplest form of the 
parametric instability model would predict resonances. 

arately the second, third, and fourth quarter of the data 
(each encompassing a period of At « 160risco) to fook 
for trends or stabihty in the PSD, while discarding the 
first quarter as potentiaUy biased by initial conditions. 
The results are summarized in Table [2l where we com- 
pare single power law models for the PSD with models 
involving a broken power law. We use the method de- 
scribed in the previous section; note that only differences 
in the log likelihood ln£ are important rather than the 
absolute magnitude of ln£. 

From this analysis, we find compelling evidence for a 
break in the power law characterizing the PSD of M. 
Compared to a single power law, the broken power fit 
is better by A In £ = 76 — 90, hence the maximum like- 
lihood ratio is at least exp(76) = lO^'^ in all three data 
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Fig. 14. — The instantaneous mass accretion rate onto the black 
hole from our canonical MHD simulation (run MHD_1). 



segments independently. The break frequency and power 
law slopes are consistent between the second and third 
data segments, but these do not match the first data seg- 
ment. The break frequencies are within a factor of two 
of, but not consistent with, the orbital frequency at the 
ISCO, logio(i'isco) = —1.79. Therefore, although there 
is a clear steepening in the power density spectrum, it 
is not possible at this point to assign a specific physical 
meaning to the break. We note that this general form of 
the PSD, i.e., an approximate power- law with curvature 
or a break at frequencies close to the ISCO orbital fre- 
quency, has been previously seen in the mass accretion 
rate of global disk simulations (Hawley & Krolik 2001, 
2002). 

We see no indications of QPOs in the M PSD. Quanti- 
tatively, we add a Lorentzian QPO to the broken power- 
law PSD model in which the QPO centroid frequency, 
full-width half-maximum, and amplitude are allowed to 
be free parameters, with the one restriction that the qual- 
ity factor of the QPO must exceed 2. We find that the 
peak power of any QPO cannot exceed 2% of the contin- 
uum power measured at the centroid of the QPO at the 
99% confidence level. 

5. DISCUSSION 
5.1. Comparison with previous numerical results 

In recent years, several groups have reported tempo- 
ral analyses of MHD disk simulations. Here, we briefly 
compare our work with some of these published results. 

Probably the most relevant previous work is that of 
Arras, Blaes & Turner (2006; hereafter ABT). These au- 
thors perform a local, shearing box MHD simulation of a 
patch of an accretion disk; this provides a controlled envi- 
ronment in which fluid modes can be characterized. ABT 
find that the MHD turbulence excites a spectrum of dis- 
tinct acoustic modes as well as radial epicyclic motions. 
However, they note a lack of distinct inertial modes (g- 
modes) and use this fact to argue against the excitation of 
trapped g- modes in global accretion disks. Our findings 
are completely in line with those of ABT, and represent 
an extension of ABT's conclusions to global simulations 
of thin accretion disks. 

There have been QPOs reported from global simula- 
tions. Kato (2004d) performed a global MHD disk sim- 
ulation in a PW potential and presented an analysis of 
quantities derived from the mass accretion rate. Through 
visual inspection of the resulting PSDs from four periods 
of the (long) simulation, he reports a pair of transient 
QPOs and a pair of QPOs that are labeled as persistent 
(although it is not clear that they are present in the PSD 
of all data segments, and the statistical significance of the 
features is unclear). The QPOs are attributed to reso- 
nances between the vertical and radial epicyclic frequen- 
cies, and it is found that these QPO pairs have frequency 
ratios that are approximately 3:2. We find no evidence 
of these resonances in our simulations. In another inter- 
esting difference, an inspection of the radially-resolved 
PSDs in Kato (2004d) reveal no signs of the local p- 
modes that seem to feature so prominently in our PSDs. 
While the reason for these discrepancies is unclear, we 
do note that the Kato (2004d) simulations have an order 
of magnitude less resolution in both the azimuthal and 
(more importantly) the vertical direction, although his 
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Fig. 15. — Comparison of the PSD and Fourier phase of the mass accretion rate without (left panel) and with (right panel) multiplication 
of the time series by an exponential in time designed to compensate for the slow loss of mass from the disk. The similarity of the two 
implies that at these frequencies the variability is dominated by intrinsic fluctuations instead of by secular changes in the disk. 



Data segment 


Model 




Parameters and uncertainties 


Maximum In C 


1 


Single PL 




F = 2.50 ±0.05 


9131 




Broken PL 


Ti 


= 0.82 ± 0.09, F2 = 2.96 ± 0.07, 
logio(l'broak) = -1.94 ±0.04 


9238 


2 


Single PL 




F = 2.70 ±0.08 


6508 




Broken PL 


Ti 


= 1.45 ± 0.11, F2 = 2.91 ± 0.10, 

logio(^broak) = -1-65 ± 0.07 


6609 


3 


Single PL 




F = 2.16 ±0.04 


9561 




Broken PL 


Ti 


= 1.47 ± 0.07, F2 = 2.89 ± 0.09, 

logio(i^break) = "1.60 ±0.04 


9630 



TABLE 2 

Model comparisons for PSD of mass accretion rate M. F is the power law index of a single power-law model for the PSD; 

P{y) tx . Fi AND F2 ARE THE TWO INDICES OBTAINED FROM A BROKEN POWER LAW; P OC V~^'^ FOR V < Z^break, P <^ FOR 

V > I'broak- All error BARS ARE ONE STANDARD DEVIATION. 



simulations do have a significantly larger computational 
domain. It is possible that the Kato (2004d) simulations 
have failed to adequately resolve the vertical dynamics 
of the thin disk. 

Chan et al. (2006) have also performed global MHD 
disk simulations in a PW potential using the pseudo- 
spectral algorithm of Chan, Psaltis, & Ozel (2005, 2006). 
Their study included a post-processing of the simulation 
to include detailed radiative transfer, and was explicitly 
targeted at understanding the variability (including the 
large amplitude flaring) of the hot accretion flow around 
the black hole at the center of the Galaxy. They found 
that the turbulence of the quiescent flow could only pro- 
duce a factor of two modulations in the observed lumi- 
nosity. To model the large amplitude flares, they in- 
troduce large density perturbations into the flow. Af- 
ter being perturbed, the disk displays a QPO with a 
frequency equal to the orbital frequency at the magne- 
tosonic point. Given that the Chan et al. study is ex- 
ploring a rather different regime of accretion than our 
present study (i.e., hot, thick accretion flows versus thin, 
cold accretion flows) , it is hard to make a direct compar- 
ison. 

Finally, SKH have performed detailed analyses of Gen- 



eral Relativistic MHD (GRMHD) simulations of disks 
performed using the code of De Villiers & Hawley (2003). 
In particular, they have studied a long (6000G'M/c ) 
simulation of a disk around a Schwarzschild black hole, 
focusing on the temporal behavior of proxy-lightcurves 
(rather than the underlying fluid properties discussed in 
this paper). It is unclear from their analysis whether 
their simulated accretion disk has excited local p-modes 
of the type that we find in this current work. SKH do 
find, though, that a proxy-lightcurve based on radiative 
transfer through the disk assuming black body emission 
and free-free absorption displays QPOs with an approxi- 
mate 3:2 ratio. However, these QPOs are transient, only 
appearing at certain times and certain viewing inclina- 
tions. 

5.2. Comments on 3:2 frequency ratios 

As described above, several authors have reported 
transient QPO pairs from MHD simulations with fre- 
quency ratio 3:2. It is tempting to interpret these as 
resonance phenomena. Here, we note that there are 
several ways that approximate 3:2 ratios can be gen- 
erated that do not necessarily involve resonances and, 
hence, one should guard against over-interpreting QPO 
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pairs. Indeed, we must remember that some sources have 
QPO frequency ratios that are inconsistent with 3:2 (e.g., 
the 67 Hz and 41 Hz QPOs from GRS 1915+105; see 
Strohmayer 2001). 

For example, consider the local vertical p-modes of 
the disk (discussed in § 12. ip . At a given radius, the 
n = 1 vertical pressure mode has a frequency of Wvcrt,! = 
•\/7 + 1 fi, where is the orbital frequency. For a gas 
pressure dominated disk in which 'y = 5/3, Wvert ,i/^ = 
1.63, and for a radiation pressure dominated disk in 
which 7 = 4/3, Wvort,!/^ = 1-53. If some unspecified 
physical process enhances emission from a given ring of 
the disk, one can imagine a situation where QPOs are 
generated at the orbital frequency and the n = 1 vertical 
pressure mode, thereby giving frequency ratios compati- 
ble with the measured ratios in several sources. Alterna- 
tively, the next lowest vertical mode that has a vertical 
velocity node in the midplane (and thus maximum vari- 
ation of pressure and density there) has a frequency of 
Wvert,3 — \/37 + 1 As a result, when 7 = 5/3 we have 
Wvert.s/wvcrt.i — 1-5, and wheu 7 — 4/3 the ratio is 1.46. 
Once again, if nature picks out a specific radius and the 
emission is modulated by the n = 1 and n = 3 modes, we 
would see QPOs with frequency ratios entirely consistent 
with the observations. 

As another example, suppose that the disk emission is 
modulated at the vertical and radial epicyclic frequen- 
cies, and that the emission is distributed in radius ac- 
cording to a standard Page & Thorne (1974) disk. The 
radial distribution of the emission then peaks close to the 
radius where the radial epicyclic frequency is a maximum 
(and hence is slowly changing with radius). One might 
then expect to see a pair of QPOs corresponding to these 
two epicyclic frequencies. The frequency ratio depends 
only weakly on the spin parameter, ranging from 1.46 at 
a/M = to 1.7 at a/M = 0.9. 

6. CONCLUSIONS 

The origin of HFQPOs remains elusive. Our simula- 
tions of geometrically-thin accretion disks have shown 
that MRI-driven MHD turbulence does not excite the 
trapped g-modes of Nowak & Wagoner (1992), even when 
those modes definitively exist in the equivalent hydrody- 
namic disk. We have also shown that MHD turbulence 
does not excite the parametric resonance instability of 
Abramowicz & Kluzniak (2001). Instead, the only dis- 
tinct modes found in our simulated MHD disks are local 
vertical p-modes and radial epicyclic oscillations. 

Clearly, the failure of all simulations to date to pro- 
duce stable QPO pairs of the type seen in GBHBs sug- 
gests that either the QPOs are too weak to be detected 
in the simulations or the models are missing some impor- 
tant physical ingredients. It is an open question whether 
the global disk modes or the parametric instabilities dis- 
cussed in this paper can be excited once one includes the 
full effects of GR close to rapidly spinning black holes 
and/or radiation physics. Indeed, the fact that HFQPOs 
are only seen in rather special spectral states (the soft in- 
termediate state; e.g., see Belloni 2006) when the accre- 
tion rate is thought to be comparable to the Eddington 
limit suggests that radiation physics, in particular, may 
well be important to the HFQPO mechanism. It is also 
noteworthy that the transition from the hard interme- 
diate state into the soft intermediate state is associated 



with powerful relativistic ejection events. Thus, another 
interesting possibility is that HFQPO production occurs 
in the black hole magnetosphere (i.e., the base of the jet) 
and not the accretion disk at all. 
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APPENDIX 

FUNDAMENTAL G-MODE FREQUENCY FOR VERY SMALL SOUND SPEEDS 

For very small sound speeds, there are simplifications that allow the frequency of the fundamental trapped g-mode in 
a hydrodynamic accretion disk to be obtained analytically. Here we base our analysis on the equations and formalism 
of NW92. 

NW92 examine the linearized equations describing the behavior of the scalar potential 

5u = bPjp, (Al) 
where is the Eulerian variation in the pressure. The radial equation for the perturbation is 



2 2 ' 



Qj.2 



{u? - 7Tr!2) {u? - K^) Si 



(A2) 



where lu is the mode angular frequency, fl is the orbital angular frequency, k is the radial epicyclic angular frequency, 
and 7 is the usual adiabatic index 7 = 5/3. T is determined by the quantization condition 

A +1/2 



(l-4B)i/2 



J + 1/2 



where A = T — ( and 



with C = (7 " 



7 LU^ 

l)/7 = 2/5. For the fundamental mode, we have j 



and eq. IA3I can be solved to obtain 



T = 0.4- 



0.2. 



(A3) 



(A4) 



(A5) 



Now we concentrate on low sound speeds, Cs ^ 1 (within this Appendix, speeds will be given in units of c; frequencies 
in units of c^/(GM); lengths in units of Vg). This implies that the mode frequency to ~ Umax, where Kmax is the 
maximum radial epicyclic frequency. In this limit, we note that the factor (w^ — jTO,"^) within eqn. IA2I is close to 



constant with radius. Let us define D = — (w — 'jTQ ) > 0, where we evaluate D by assuming uj — Kmax- 
We also note that near the maximum, the radial epicyclic frequency has a parabolic form, — K,^„ax ~ E{r- 



where r^ax is the radius at which k ~ k„ 



Qj.2 



The differential equation then becomes 
D , ^ ,9 X D E 



('^max/'^ - 1) + — — 

ci ci LW^ 



5u 



(A6) 
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This has the form of a harmonic osciUator, so we try a solution of the type 

Su cx e-^^('-'-— 



meaning that 

This yields the conditions 



e 2 



,2 

's 

D E 



Defining x = these two conditions can be combined (eliminating C) to yield 



Solving for x gives 



(A7) 
(A8) 

(A9) 

(AlO) 
(All) 

(A12) 
(A13) 



Since uj < Kmax, we choose the positive sign. For small Cs, the first term in the square root dominates, and the lowest 
order in Cc gives 



1 



1 



This implies finally that 



-VE/D 



(A14) 



(A15) 



This clearly demonstrates that the fractional difference of the mode frequency from the radial epicyclic maximum is 
linear in c^. 

Evaluating this expression for the PW potential, we obtain k^^x = 1-202 x lO^-^ (at r = Tmax = 4 + 2V3 = 7.464), 
and 



= 2.269cs . 

^max 

In contrast, the Nowak & Wagoner pseudo-Newtonian potential, 

1 3 12 

4'nw = 

yields k^^^^ = 8.607 x lO"'' (at r^ax = 7.746) and 



5.621cs 



(A16) 



(A17) 



(A18) 



For both potentials, the analytic frequencies agree extremely well with direct numerical solutions to eqn. IA2I 



